CMIP6 prep
Camille
2024-02-19
library(raster)
library(plyr)
library(dplyr)
library(ggplot2)
library(viridisLite)
library(colorspace)
library(RColorBrewer)
library(sf)
library(stringr)
library(ggpubr)
library(kableExtra)
memory.limit(size = 160000)
## [1] Inf
display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")1 Data
land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land
elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif")
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))
ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326"))
ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))tas.all_annual <- read.table("CMIP6/tas/tas.all_annual.txt")
pr.all_annual <- read.table("CMIP6/pr/pr.all_annual.txt")
sfcWind.all_annual <- read.table("CMIP6/sfcWind/sfcWind.all_annual.txt")
hfls.all_annual <- read.table("CMIP6/hfls/hfls.all_annual.txt")
hfss.all_annual <- read.table("CMIP6/hfss/hfss.all_annual.txt")
hurs.all_annual <- read.table("CMIP6/hurs/hurs.all_annual.txt")
cmip6_annual <- tas.all_annual %>% merge(pr.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_annual, "cmip6_annual.txt")tas.all_january <- read.table("CMIP6/tas/tas.all_january.txt")
pr.all_january <- read.table("CMIP6/pr/pr.all_january.txt")
sfcWind.all_january <- read.table("CMIP6/sfcWind/sfcWind.all_january.txt")
hfls.all_january <- read.table("CMIP6/hfls/hfls.all_january.txt")
hfss.all_january <- read.table("CMIP6/hfss/hfss.all_january.txt")
hurs.all_january <- read.table("CMIP6/hurs/hurs.all_january.txt")
cmip6_january <- tas.all_january %>% merge(pr.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_january, "cmip6_january.txt")tas.all_february <- read.table("CMIP6/tas/tas.all_february.txt")
pr.all_february <- read.table("CMIP6/pr/pr.all_february.txt")
sfcWind.all_february <- read.table("CMIP6/sfcWind/sfcWind.all_february.txt")
hfls.all_february <- read.table("CMIP6/hfls/hfls.all_february.txt")
hfss.all_february <- read.table("CMIP6/hfss/hfss.all_february.txt")
hurs.all_february <- read.table("CMIP6/hurs/hurs.all_february.txt")
cmip6_february <- tas.all_february %>% merge(pr.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_february, "cmip6_february.txt")tas.all_march <- read.table("CMIP6/tas/tas.all_march.txt")
pr.all_march <- read.table("CMIP6/pr/pr.all_march.txt")
sfcWind.all_march <- read.table("CMIP6/sfcWind/sfcWind.all_march.txt")
hfls.all_march <- read.table("CMIP6/hfls/hfls.all_march.txt")
hfss.all_march <- read.table("CMIP6/hfss/hfss.all_march.txt")
hurs.all_march <- read.table("CMIP6/hurs/hurs.all_march.txt")
cmip6_march <- tas.all_march %>% merge(pr.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_march, "cmip6_march.txt")tas.all_april <- read.table("CMIP6/tas/tas.all_april.txt")
pr.all_april <- read.table("CMIP6/pr/pr.all_april.txt")
sfcWind.all_april <- read.table("CMIP6/sfcWind/sfcWind.all_april.txt")
hfls.all_april <- read.table("CMIP6/hfls/hfls.all_april.txt")
hfss.all_april <- read.table("CMIP6/hfss/hfss.all_april.txt")
hurs.all_april <- read.table("CMIP6/hurs/hurs.all_april.txt")
cmip6_april <- tas.all_april %>% merge(pr.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_april, "cmip6_april.txt")tas.all_may <- read.table("CMIP6/tas/tas.all_may.txt")
pr.all_may <- read.table("CMIP6/pr/pr.all_may.txt")
sfcWind.all_may <- read.table("CMIP6/sfcWind/sfcWind.all_may.txt")
hfls.all_may <- read.table("CMIP6/hfls/hfls.all_may.txt")
hfss.all_may <- read.table("CMIP6/hfss/hfss.all_may.txt")
hurs.all_may <- read.table("CMIP6/hurs/hurs.all_may.txt")
cmip6_may <- tas.all_may %>% merge(pr.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_may, "cmip6_may.txt")tas.all_june <- read.table("CMIP6/tas/tas.all_june.txt")
pr.all_june <- read.table("CMIP6/pr/pr.all_june.txt")
sfcWind.all_june <- read.table("CMIP6/sfcWind/sfcWind.all_june.txt")
hfls.all_june <- read.table("CMIP6/hfls/hfls.all_june.txt")
hfss.all_june <- read.table("CMIP6/hfss/hfss.all_june.txt")
hurs.all_june <- read.table("CMIP6/hurs/hurs.all_june.txt")
cmip6_june <- tas.all_june %>% merge(pr.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_june, "cmip6_june.txt")tas.all_july <- read.table("CMIP6/tas/tas.all_july.txt")
pr.all_july <- read.table("CMIP6/pr/pr.all_july.txt")
sfcWind.all_july <- read.table("CMIP6/sfcWind/sfcWind.all_july.txt")
hfls.all_july <- read.table("CMIP6/hfls/hfls.all_july.txt")
hfss.all_july <- read.table("CMIP6/hfss/hfss.all_july.txt")
hurs.all_july <- read.table("CMIP6/hurs/hurs.all_july.txt")
cmip6_july <- tas.all_july %>% merge(pr.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_july, "cmip6_july.txt")tas.all_august <- read.table("CMIP6/tas/tas.all_august.txt")
pr.all_august <- read.table("CMIP6/pr/pr.all_august.txt")
sfcWind.all_august <- read.table("CMIP6/sfcWind/sfcWind.all_august.txt")
hfls.all_august <- read.table("CMIP6/hfls/hfls.all_august.txt")
hfss.all_august <- read.table("CMIP6/hfss/hfss.all_august.txt")
hurs.all_august <- read.table("CMIP6/hurs/hurs.all_august.txt")
cmip6_august <- tas.all_august %>% merge(pr.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_august, "cmip6_august.txt")tas.all_september <- read.table("CMIP6/tas/tas.all_september.txt")
pr.all_september <- read.table("CMIP6/pr/pr.all_september.txt")
sfcWind.all_september <- read.table("CMIP6/sfcWind/sfcWind.all_september.txt")
hfls.all_september <- read.table("CMIP6/hfls/hfls.all_september.txt")
hfss.all_september <- read.table("CMIP6/hfss/hfss.all_september.txt")
hurs.all_september <- read.table("CMIP6/hurs/hurs.all_september.txt")
cmip6_september <- tas.all_september %>% merge(pr.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_september, "cmip6_september.txt")tas.all_october <- read.table("CMIP6/tas/tas.all_october.txt")
pr.all_october <- read.table("CMIP6/pr/pr.all_october.txt")
sfcWind.all_october <- read.table("CMIP6/sfcWind/sfcWind.all_october.txt")
hfls.all_october <- read.table("CMIP6/hfls/hfls.all_october.txt")
hfss.all_october <- read.table("CMIP6/hfss/hfss.all_october.txt")
hurs.all_october <- read.table("CMIP6/hurs/hurs.all_october.txt")
cmip6_october <- tas.all_october %>% merge(pr.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_october, "cmip6_october.txt")tas.all_november <- read.table("CMIP6/tas/tas.all_november.txt")
pr.all_november <- read.table("CMIP6/pr/pr.all_november.txt")
sfcWind.all_november <- read.table("CMIP6/sfcWind/sfcWind.all_november.txt")
hfls.all_november <- read.table("CMIP6/hfls/hfls.all_november.txt")
hfss.all_november <- read.table("CMIP6/hfss/hfss.all_november.txt")
hurs.all_november <- read.table("CMIP6/hurs/hurs.all_november.txt")
cmip6_november <- tas.all_november %>% merge(pr.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_november, "cmip6_november.txt")tas.all_december <- read.table("CMIP6/tas/tas.all_december.txt")
pr.all_december <- read.table("CMIP6/pr/pr.all_december.txt")
sfcWind.all_december <- read.table("CMIP6/sfcWind/sfcWind.all_december.txt")
hfls.all_december <- read.table("CMIP6/hfls/hfls.all_december.txt")
hfss.all_december <- read.table("CMIP6/hfss/hfss.all_december.txt")
hurs.all_december <- read.table("CMIP6/hurs/hurs.all_december.txt")
cmip6_december <- tas.all_december %>% merge(pr.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(sfcWind.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfls.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hfss.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(hurs.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
merge(elev.df, by = c("lon", "lat"))
write.table(cmip6_december, "cmip6_december.txt")2 Aridity Index
cmip6 <- read.table("cmip6_annual.txt")
breaks.martonne <- c(0,5,10,20,"30","40",Inf)
cat.martonne <- c("[0,5]" = "Desert", "(5,10]"= "Arid", "(10,20]" = "Semi-arid","(20,30]" = "Temperate", "(30,40]" = "Humid", "(40,Inf]" = "Forest")
breaks.unesco <- c(-Inf,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(-Inf,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid")
col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")2.1 De Martonne
cmip6$AIm <- with(cmip6,pr*60*60*24*365/(tas-273.15+10)) # pr in mm/y, tas in Ceslsius
cmip6$cat.AIm <- cmip6$AIm %>% cut(breaks.martonne, include.lowest = T) %>% revalue(cat.martonne)
ggplot(subset(cmip6, period == "1970_2000")) + geom_raster(aes(x=lon, y = lat, fill = cat.AIm))+
scale_fill_manual(values = col.martonne)+
theme_void()+
theme(legend.position = "bottom")
2.2 FAO
ea from Rh:
ea = 6.108(RH/100)exp(17.27*T/(T+237.3)) Refs https://www.nature.com/articles/s41598-023-40499-6 : Allen and Tetens
2.2.1 ET0 by month
id_vars <- c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")
cmip6_annual <- mutate(read.table("cmip6_annual.txt"),
t = tas - 273.15,
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_annual = ((0.408*delta*Rn+gamma*900/((t)+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_01 <- mutate(read.table("cmip6_january.txt"),
t = tas - 273.15,
pr_01 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_01 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) # Evapotranspiration, G considered to be negligible. in mm/day. Wind speed at 10 m is converted to wind speed at 2 m by multiplying by 0.748
cmip6_02 <- mutate(read.table("cmip6_february.txt"),
t = tas - 273.15,
pr_02 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_02 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_03 <- mutate(read.table("cmip6_march.txt"),
t = tas - 273.15,
pr_03 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_03 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_04 <- mutate(read.table("cmip6_april.txt"),
t = tas - 273.15,
pr_04 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_04 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_05 <- mutate(read.table("cmip6_may.txt"),
t = tas - 273.15,
pr_05 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_05 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_06 <- mutate(read.table("cmip6_june.txt"),
t = tas - 273.15,
pr_06 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_06 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_07 <- mutate(read.table("cmip6_july.txt"),
t = tas - 273.15,
pr_07 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_07 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_08 <- mutate(read.table("cmip6_august.txt"),
t = tas - 273.15,
pr_08 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_08 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_09 <- mutate(read.table("cmip6_september.txt"),
t = tas - 273.15,
pr_09 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_09 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_10 <- mutate(read.table("cmip6_october.txt"),
t = tas - 273.15,
pr_10 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_10 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_11 <- mutate(read.table("cmip6_november.txt"),
t = tas - 273.15,
pr_11 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_11 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
cmip6_12 <- mutate(read.table("cmip6_december.txt"),
t = tas - 273.15,
pr_12 = pr*60*60*24, # mm/day
Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
P = 101.3 * ((293-0.0065*z)/293)^(5.26),
es = 0.61078*exp(17.27*t/(t+237.3)),
ea = (hurs/100)*es,
delta = 4098*es/(t+237.3)^2, #slope vapor pressure
gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
ET0_12 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))
et0_monthly <- select(cmip6_annual, id_vars) %>%
merge(select(cmip6_01, c(id_vars, "ET0_01", "pr_01")), by = id_vars) %>%
merge(select(cmip6_02, c(id_vars, "ET0_02", "pr_02")), by = id_vars) %>%
merge(select(cmip6_03, c(id_vars, "ET0_03", "pr_03")), by = id_vars) %>%
merge(select(cmip6_04, c(id_vars, "ET0_04", "pr_04")), by = id_vars) %>%
merge(select(cmip6_05, c(id_vars, "ET0_05", "pr_05")), by = id_vars) %>%
merge(select(cmip6_06, c(id_vars, "ET0_06", "pr_06")), by = id_vars) %>%
merge(select(cmip6_07, c(id_vars, "ET0_07", "pr_07")), by = id_vars) %>%
merge(select(cmip6_08, c(id_vars, "ET0_08", "pr_08")), by = id_vars) %>%
merge(select(cmip6_09, c(id_vars, "ET0_09", "pr_09")), by = id_vars) %>%
merge(select(cmip6_10, c(id_vars, "ET0_10", "pr_10")), by = id_vars) %>%
merge(select(cmip6_11, c(id_vars, "ET0_11", "pr_11")), by = id_vars) %>%
merge(select(cmip6_12, c(id_vars, "ET0_12", "pr_12")), by = id_vars)
et0_monthly$spr <- with(et0_monthly, pr_01*31 + pr_02 * 28 + pr_03 * 31 + pr_04 * 30 + pr_05 * 31 + pr_06 * 30 +
pr_07 * 31 + pr_08 * 31 + pr_09 * 30 + pr_10 * 31 + pr_11 * 30 + pr_12 * 31)
et0_monthly$sET0 <- with(et0_monthly, ET0_01*31 + ET0_02 * 28 + ET0_03 * 31 + ET0_04 * 30 + ET0_05 * 31 + ET0_06 * 30 +
ET0_07 * 31 + ET0_08 * 31 + ET0_09 * 30 + ET0_10 * 31 + ET0_11 * 30 + ET0_12 * 31)
cmip6 <- merge(cmip6_annual, et0_monthly, by = id_vars)
cmip6$AI <- with(cmip6, abs(spr/sET0))
cmip6$AI_annual <- with(cmip6, abs((pr*60*60*24*365)/(ET0_annual*365)))
breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6$cat.AI[which(cmip6$sET0 < 400)] <- "Cold"
cmip6$cat.AI %>% unique()
cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>%
revalue(cat.unesco)
cmip6$cat.AI_annual <- cmip6$AI_annual %>% cut(breaks = breaks.unesco) %>%
revalue(cat.unesco)
cmip6$cat.AI_annual[which(cmip6$ET0_annual*365 < 400)] <- "Cold"
write.table(cmip6, "cmip6_AI.txt")
cmip6 <- read.table("cmip6_AI.txt")
cmip6_AI <- cmip6 %>% subset(period == "1970_2000") %>%
group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>%
dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(pr, na.rm = T) * 60 * 60 * 24 * 365,
sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>%
ungroup() %>%
mutate(source = "CMIP6")%>%
select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))
write.table(cmip6_AI, "cmip6ref.txt")3 Maps ETO by model
cmip6 <- read.table("cmip6_AI.txt")3.1 1850-1880
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "1850_1880" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
3.2 1970-2000
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
library(ggrepel)
pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
select(c("source", "cat.AI")) %>%
reshape2::melt(id.vars = c("source")) %>%
group_by(source, value) %>%
summarise(count = n()) %>%
mutate(perc = round(count/sum(count)*100, 1),
lab.y = rev(cumsum(rev(count))) - count*0.5) %>%
ungroup()
pie_list <- list()
for(i in unique(pie$source)){
p <- ggplot(subset(pie, source == i))+
geom_bar(aes(x = "", y = count, fill = value), width = 1, stat = "identity")+
coord_polar("y", start = 0)+
scale_fill_manual(values = col.cat)+
geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
labs(fill = "", title = i)+
theme_void()+theme(legend.position = "right")
pie_list[[i]] <- p
}
ggpubr::ggarrange(plotlist = pie_list, nrow = 4, ncol = 4, common.legend = T, legend = "bottom")
multimodel <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(cat.AI) %>%
summarise(count = n()) %>%
mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
tidyr::pivot_wider(names_from = cat.AI, values_from = percent) %>%
mutate(source = "Multimodel average","Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
tab <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(source, cat.AI) %>%
summarise(count = n()) %>%
ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
reshape2::dcast(source~cat.AI) %>%
mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>%
rbind(multimodel)
knitr::kable(tab, digits = 2, escape = F) %>%
kable_styling(bootstrap_options = "bordered") %>%
column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(14, bold = T) | source | Hyperarid | Arid | Semi-arid | Dry subhumid | Sum drylands | Humid | Cold | NA |
|---|---|---|---|---|---|---|---|---|
| CAS-ESM2 | 3.6 | 6.0 | 9.2 | 4.4 | 23.2 | 26.0 | 48.2 | 2.6 |
| CESM | 3.7 | 7.4 | 9.8 | 3.6 | 24.5 | 24.3 | 48.5 | 2.6 |
| CMCC | 3.6 | 5.2 | 7.0 | 3.8 | 19.6 | 29.4 | 48.5 | 2.6 |
| CMCC-ESM2 | 4.5 | 5.5 | 5.3 | 1.8 | 17.1 | 17.0 | 63.3 | 2.6 |
| CNRM | 4.0 | 8.5 | 7.9 | 4.5 | 24.9 | 26.6 | 45.9 | 2.6 |
| EC-Earth3 | 5.9 | 6.6 | 7.5 | 3.1 | 23.1 | 24.8 | 49.6 | 2.6 |
| FGOALS | 4.1 | 7.6 | 7.0 | 3.4 | 22.1 | 17.5 | 57.9 | 2.6 |
| GFDL-ESM4 | 4.0 | 6.2 | 6.2 | 2.8 | 19.2 | 24.6 | 53.6 | 2.6 |
| INM | 3.1 | 5.9 | 9.2 | 4.0 | 22.2 | 27.9 | 47.4 | 2.6 |
| INM-CM5 | 2.3 | 5.9 | 8.4 | 3.9 | 20.5 | 29.9 | 47.0 | 2.6 |
| MPI | 6.0 | 7.7 | 7.1 | 3.0 | 23.8 | 23.7 | 49.9 | 2.6 |
| MRI | 4.3 | 7.5 | 6.4 | 2.7 | 20.9 | 27.0 | 49.5 | 2.6 |
| NorESM-2-MM | 3.6 | 7.0 | 11.9 | 3.9 | 26.4 | 23.3 | 47.7 | 2.6 |
| Multimodel average | 4.0 | 6.7 | 7.9 | 3.5 | 22.1 | 24.8 | 50.5 | 2.6 |
3.3 1985-2010
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "1985_2015" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "1985-2010", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
3.4 SSP 2-4.5
3.4.1 2030-2060
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
3.4.2 2070-2100
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
3.5 SSP 3-7.0
3.5.1 2030-2060
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
3.5.2 2070-2100
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
3.6 SSP 5-8.5
3.6.1 2030-2060
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
3.6.2 2070-2100
map_list <- list()
for(i in unique(cmip6$source)){
g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
theme_void()
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")
4 Agreement between models
4.1 Summary dataframe
cmip6s <- cmip6 %>%
group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, period, z) %>%
dplyr::summarise(t.mean = mean(t, na.rm = T), t.sd = sd(t, na.rm = T), # t in °C
pr.mean = mean(spr, na.rm = T) , pr.sd = sd(spr, na.rm = T), # pr in mm/y
sfcWind.mean = mean(sfcWind, na.rm = T), sfcWind.sd = sd(sfcWind, na.rm = T),
Rn.mean = mean(Rn, na.rm = T), Rn.sd = sd(Rn, na.rm = T), # in MJ/m2/day
ea.mean = mean(ea, na.rm = T), ea.sd = sd(ea, na.rm = T), # in kPa
ET0 = mean(sET0, na.rm = T), ET0.sd = sd(sET0, na.rm = T), # in mm/y
AI.mean = mean(AI, na.rm = T), AI.sd = sd(AI, na.rm = T),
hyperarid = sum(cat.AI == "Hyperarid"),
arid = sum(cat.AI == "Arid"),
semiarid = sum(cat.AI == "Semi-arid"),
drysubhumid = sum(cat.AI == "Dry subhumid") ,
humid = sum(cat.AI == "Humid"),
cold = sum(cat.AI == "Cold"),
na = sum(is.na(cat.AI)),
nb.models = n(),
nb.maj = max(hyperarid,arid,semiarid,drysubhumid,humid,cold,na),
prop.maj = nb.maj/nb.models,
) %>% ungroup() %>% as.data.frame()
cmip6s$cat.maj <- colnames(select(cmip6s, c("hyperarid","arid",'semiarid',"drysubhumid",'humid',"cold","na")))[max.col(select(cmip6s, c("hyperarid","arid",'semiarid',"drysubhumid",'humid',"cold","na")),ties.method = "first")]
write.table(cmip6s, "cmip6s.txt")cmip6s <- read.table("cmip6s.txt")
breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
cmip6s$cat.AI <- cmip6s$AI.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI[which(cmip6s$ET0 < 400)] <- "Cold"
# getAIref <- function(i) {
# AIrefi = subset(cmip6s, lat == cmip6s$lat[i] & lon == cmip6s$lon[i] & period =="1970_2000") %>% pull(AI.mean)
# return(AIrefi)
# }
#
# AI.ref <- sapply(as.numeric(rownames(cmip6s)), FUN = getAIref)
# saveRDS(AI.ref, "AI.ref.rds")
cmip6s$AI.ref <- readRDS("AI.ref.rds")
# getcatAIref <- function(i) {
# catAIrefi = subset(cmip6s, lat == cmip6s$lat[i] & lon == cmip6s$lon[i] & period =="1970_2000") %>% pull(cat.AI)
# return(catAIrefi)
# }
#
# catAIref.vect <- sapply(as.numeric(rownames(cmip6s)), FUN = getcatAIref)
# saveRDS(catAIref, "cat.AI.ref.rds")
cmip6s$catAIref <- readRDS("catAIref.vect.rds")
# getcatmajref <- function(i) {
# catmajrefi = subset(cmip6s, lat == cmip6s$lat[i] & lon == cmip6s$lon[i] & period =="1970_2000") %>% pull(cat.maj)
# return(catmajrefi)
# }
#
# cmip6.catmajref <- sapply(as.numeric(rownames(cmip6s)), FUN = getcatmajref)
# saveRDS(cmip6.catmajref, "cmip6.catmajref.rds")
cmip6s$catmajref <- readRDS("cmip6.catmajref.rds")
cmip6s$anomalies <- with(cmip6s, AI.mean - AI.ref)
write.table(cmip6s, "cmip6s2.txt")cmip6s <- read.table("cmip6s2.txt")
# attributes a change from the studied period compared to the reference period 1970-2000. Read as: grid cell x is humid in 2030-2060 and Arid in 1970-2000, so it is WETTER
#cmip6s$change <- with(cmip6s, paste(catAIref, cat.AI, collapse = " ", sep = " to "))
cmip6s$dryer <- with(cmip6s, ifelse(catAIref == "Arid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter",
ifelse(catAIref == "Arid" & cat.AI == "Hyperarid", "dryer",
ifelse(catAIref == "Semi-arid" & cat.AI %in% c("Cold","Humid","Dry subhumid"), "wetter",
ifelse(catAIref == "Semi-arid" & cat.AI %in% c("Hyperarid","Arid"), "dryer",
ifelse(catAIref == "Dry subhumid" & cat.AI %in% c("Cold","Humid"), "wetter",
ifelse(catAIref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
ifelse(catAIref == "Humid" & cat.AI == "Cold","wetter",
ifelse(catAIref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
ifelse(catAIref == "Cold" & cat.AI %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
ifelse(catAIref == "Hyperarid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
NA)))))))))))
cmip6s$catmajref <- revalue(cmip6s$catmajref, c("hyperarid" = "Hyperarid", "arid" = "Arid", "semiarid" = "Semi-arid", "drysubhumid" = "Dry subhumid", "humid" = "Humid", "cold" = "Cold"))
cmip6s$cat.maj <- revalue(cmip6s$cat.maj, c("hyperarid" = "Hyperarid", "arid" = "Arid", "semiarid" = "Semi-arid", "drysubhumid" = "Dry subhumid", "humid" = "Humid", "cold" = "Cold"))
cmip6s$dryer.maj <- with(cmip6s, ifelse(catmajref == "Arid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter",
ifelse(catmajref == "Arid" & cat.maj == "Hyperarid", "dryer",
ifelse(catmajref == "Semi-arid" & cat.maj %in% c("Cold","Humid","Dry subhumid"), "wetter",
ifelse(catmajref == "Semi-arid" & cat.maj %in% c("Hyperarid","Arid"), "dryer",
ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Cold","Humid"), "wetter",
ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
ifelse(catmajref == "Humid" & cat.maj == "Cold","wetter",
ifelse(catmajref == "Humid" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
ifelse(catmajref == "Cold" & cat.maj %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
ifelse(catmajref == "Hyperarid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
NA)))))))))))
4.2 Compare cat mean and cat maj
4.2.1 Where are cat mean and cat maj different
cmip6s$same.catAI <- with(cmip6s, ifelse(cat.maj == cat.AI, "same", "diff"))
cmip6s %>% group_by(model, period) %>% summarise(maj.equal.mean = table(same.catAI)[2], maj.diff.mean = table(same.catAI)[1], prop.same = round(maj.equal.mean/n()*100,1))
## # A tibble: 9 × 5
## # Groups: model [4]
## model period maj.equal.mean maj.diff.mean prop.same
## <chr> <chr> <int> <int> <dbl>
## 1 SSP245 2030_2060 19884 1798 89.3
## 2 SSP245 2070_2100 19890 1792 89.3
## 3 SSP370 2030_2060 19988 1694 89.8
## 4 SSP370 2070_2100 19911 1771 89.4
## 5 SSP585 2030_2060 19799 1883 88.9
## 6 SSP585 2070_2100 19796 1886 88.9
## 7 historical 1850_1880 19907 1775 89.4
## 8 historical 1970_2000 19880 1802 89.3
## 9 historical 1985_2015 19936 1746 89.6
map_list <- list()
for(i in unique(cmip6s$model)){
df <- subset(cmip6s, model == i)
for(j in unique(df$period)){
index <- paste(i, j, sep = "")
g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.AI))+
borders()+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, j , sep = " "), fill = "")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[index]] <- g
}
}
ggarrange(plotlist = map_list, nrow = 3, ncol = 3, common.legend = T)
map_list2 <- list()
for(i in unique(cmip6s$model)){
df <- subset(cmip6s, model == i)
for(j in unique(df$period)){
index <- paste(i, j, sep = "")
g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.maj))+
borders()+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste(i, j , sep = " "), fill = "")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list2[[index]] <- g
}
}
ggarrange(plotlist = map_list2, nrow = 3, ncol = 3, common.legend = T)
cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))
cmip6s$cat.maj <- factor(cmip6s$cat.maj, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))
ggplot(subset(cmip6s, cat.AI != "Cold"))+
geom_col(aes(x= period, y = cat.AI, fill = cat.AI), just = -0.2, width = 0.3)+
geom_col(aes(x = period, y = cat.maj, fill = cat.maj), just = 1.2, width = 0.3)+
scale_fill_manual(values = col.cat)+
labs(fill = "", y ="")+
facet_grid(rows = vars(model), switch = "y")+
theme_minimal()+theme(axis.text.y = element_blank())
tab <- rbind(table(cmip6s$cat.AI)/length(cmip6s$cat.AI)*100, table(cmip6s$cat.maj)/length(cmip6s$cat.maj)*100) %>% t() %>% as.data.frame() %>% setNames(c("Cat mean", "Cat maj"))
knitr::kable(tab, digits =2)| Cat mean | Cat maj | |
|---|---|---|
| Hyperarid | 3.51 | 4.18 |
| Arid | 6.25 | 7.43 |
| Semi-arid | 7.63 | 9.19 |
| Dry subhumid | 3.39 | 0.84 |
| Humid | 28.75 | 28.53 |
| Cold | 47.90 | 47.23 |
4.2.2 Plot cat maj and cat mean for reference period
ggpubr::ggarrange(plotlist = list(
ggplot() +
geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat, fill = cat.maj))+
geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 1, col = "grey20")+
scale_fill_manual(values = col.cat)+
labs(title = "Category chosen in majority by models", fill = "")+
ylim(-55,90)+
theme_void(base_size = 15)+
theme(legend.position = "none"),
ggplot() +
geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat, fill = cat.AI))+
geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 1, col = "grey20")+
scale_fill_manual(values = col.cat)+
labs(title = "Multimodel average categories of AI", fill = "")+
ylim(-55,90)+
theme_void(base_size = 15)+
theme(legend.position = "none")),
ncol = 2)
map_list <- list()
for(i in c("1850_1880","1970_2000","1985_2015")){
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat, fill = cat.maj))+
geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), shape = "'", col = "black")+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste("Aridity category that is predicted\nby most models in",i , sep = " "), fill = "")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[i]] <- g
}
map_list2 <- list()
for(i in c("1850_1880","1970_2000","1985_2015")){
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat, fill = cat.AI))+
geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), shape = "'", col = "black")+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(title = paste("Aridity category computed\nfor mean of models in",i , sep = " "), fill = "")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list2[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, common.legend = T)
ggpubr::ggarrange(plotlist = map_list2, ncol = 3, common.legend = T)
map_list <- list()
for(i in c("2030_2060","2070_2100")){
for (j in c("SSP245","SSP370","SSP585")){
index <- paste(i, j, sep = " , ")
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == j), aes(x=lon, y = lat, fill = cat.maj))+
# geom_point(data = subset(cmip6s, period == i & model == "historical" & change.cat == 1 & nb.maj > 10), aes(x= lon, y = lat), shape = "'")+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(fill = "Cat AI maj", title = paste(j, i, sep = ","))+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[index]] <- g
}}
map_list2 <- list()
for(i in c("2030_2060","2070_2100")){
for (j in c("SSP245","SSP370","SSP585")){
index <- paste(i, j, sep = " , ")
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == j), aes(x=lon, y = lat, fill = cat.AI))+
# geom_point(data = subset(cmip6s, period == i & model == "historical" & change.cat == 1 & nb.maj > 10), aes(x= lon, y = lat), shape = "'")+
scale_fill_manual(values = col.cat, na.translate = F)+
labs(fill = "Cat AI mean", title = paste(j, i, sep = ","))+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list2[[index]] <- g
}}
ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 2, common.legend = T)
ggpubr::ggarrange(plotlist = map_list2, ncol = 3, nrow = 2, common.legend = T)
5 Anomalies
v<-NULL
v2<-NULL
for (N in 5:30) {
v<-c(v,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N*2)))# proba d'avoir 2/3 de même signe à partir de N modèles, lorsque N change
v2<-c(v2,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N)))# proba d'avoir 2/3 de signe >0 à partir de N modèles, lorsque N change
}
plot(5:30,v,type="l",ylab="Proba(>2/3 models...)",xlab="Nb of models")
lines(5:30,v2,col="red")
5.1 Differences in AI
ggpubr::ggarrange(plotlist = list(
ggplot(subset(cmip6s, anomalies < 1 & anomalies > -1 & model == "historical"))+
geom_boxplot(aes(x=period, y = - anomalies, col = period), outlier.size = 1)+
scale_color_manual(values = c("#1A9850","#66BD63","#A6D96A"))+
labs(x = "", y = "difference of AI between each period and\nthe reference period 1970-2000")+
facet_grid(cols = vars(model))+
theme_minimal() + theme(legend.position = "none"),
ggplot(subset(cmip6s, anomalies < 1 & anomalies > -1 & model != "historical"))+
geom_boxplot(aes(x=period, y = anomalies, col = period), outlier.size = 1)+
scale_color_manual(values = c("#FDAE61","#F46D43"))+
scale_y_continuous(position = "right")+
labs(y = "", x = "")+
facet_grid(cols = vars(model))+
theme_minimal() + theme(legend.position = "none")),
widths = c(1.5,3),
ncol = 2, nrow = 1
)
5.1.1 Strictly wetter or dryer (reference period: 1970-2000)
map_list <- list()
for(i in c("1850_1880","1985_2015")){
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical" & anomalies > 0), aes(x=lon, y = lat, fill = "wetter"))+
geom_raster(data = subset(cmip6s, period == i & model == "historical" & anomalies < 0), aes(x=lon, y = lat, fill = "dryer"))+
borders()+
scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
labs(title = paste("AI changes between",i , "and 1970-2000", sep = " "), fill = "")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[i]] <- g
}
ggarrange(plotlist = map_list, ncol = 2)
map_list <- list()
for(i in c("2030_2060","2070_2100")){
for (j in c("SSP245","SSP370","SSP585")){
index <- paste(i, j, sep = " , ")
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == j & anomalies > 0), aes(x=lon, y = lat, fill = "wetter"))+
geom_raster(data = subset(cmip6s, period == i & model == j & anomalies < 0), aes(x=lon, y = lat, fill = "dryer"))+
borders()+
scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
labs(fill = "Compared to 1970-2000")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[index]] <- g
}}
ggpubr::ggarrange(plotlist = map_list,ncol = 3, nrow = 2, common.legend = T, legend = "bottom", labels = c("SSP 2-4.5", "SSP 3-7.0", "SSP 5.85"))
5.1.2 Wetter or dryer by class (reference period: 1970-2000)
ano.quantiles <- cmip6s %>% group_by(model, period) %>% summarise(q10 = quantile(anomalies, probs = 0.1, na.rm =T),
q25 = quantile(anomalies, probs = 0.25, na.rm =T),
q50 = quantile(anomalies, probs = 0.5, na.rm =T),
q75 = quantile(anomalies, probs = 0.75, na.rm =T),
q90 = quantile(anomalies, probs = 0.9, na.rm =T))
ano.quantiles <- cmip6s %>% summarise(q10 = quantile(anomalies, probs = 0.1, na.rm =T),
q25 = quantile(anomalies, probs = 0.25, na.rm =T),
q50 = quantile(anomalies, probs = 0.5, na.rm =T),
q75 = quantile(anomalies, probs = 0.75, na.rm =T),
q90 = quantile(anomalies, probs = 0.9, na.rm =T))
qbreaks <- c(-5, -1,-0.1, -0.01, 0, 0.01, 0.1, 1,5)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")“Anomalies” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.
AS OF NOW THE SAME QUANTILES ARE USED FOR ALL MAPS.
map_list <- list()
for(i in c("1850_1880","1985_2015")){
#quant = subset(ano.quantiles, model == "historical" & period == i)
#b = quant[1,c(3:7)] %>% as.numeric
b = qbreaks
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat, fill = anomalies))+
borders()+
binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
guide = guide_legend(label.theme = element_text(angle = 0)))+
labs(title = i, fill = "")+
theme_void()+ylim(-55,90)+
theme(legend.position = "right")
map_list[[i]] <- g
}
ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")
map_list <- list()
for(i in c("2030_2060","2070_2100")){
for (j in c("SSP245","SSP370","SSP585")){
index <- paste(i, j, sep = " , ")
b = qbreaks
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == j), aes(x=lon, y = lat, fill = anomalies))+
borders()+
binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
guide = guide_legend(label.theme = element_text(angle = 0)))+
labs(fill = "Compared to 1970-2000")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[index]] <- g
}}
ggpubr::ggarrange(plotlist = map_list,ncol = 3, nrow = 2, common.legend = T, legend = "bottom", labels = c("SSP 2-4.5", "SSP 3-7.0", "SSP 5.85"))
5.2 Difference between future and reference category of AI (using cat maj)
5.2.1 Dryer or wetter category
cmip6s$change.catmaj <- with(cmip6s, ifelse(cat.maj == catmajref, 0, 1))
map_list <- list()
for(i in c("1850_1880","1985_2015")){
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1), aes(x=lon, y = lat, fill = dryer.maj))+
borders()+
scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
labs(title = paste("AI changes between",i , "and 1970-2000", sep = " "), fill = "")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[i]] <- g
}
ggarrange(plotlist = map_list, ncol = 2, common.legend = T)
map_list <- list()
for(i in c("2030_2060","2070_2100")){
for (j in c("SSP245","SSP370","SSP585")){
index <- paste(i, j, sep = " , ")
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == j), aes(x=lon, y = lat, fill = dryer.maj))+
geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & nb.maj > 10), aes(x= lon, y = lat, fill = "'"))+
borders()+
scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
labs(fill = "Compared to 1970-2000")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom")
map_list[[index]] <- g
}}
ggpubr::ggarrange(plotlist = map_list,ncol = 3, nrow = 2, common.legend = T, legend = "bottom", labels = c("SSP 2-4.5", "SSP 3-7.0", "SSP 5.85"))
### All differences
5.2.2 Detailed changes and colors
cmip6s$change.catmaj <- with(cmip6s, ifelse(cat.maj == catmajref, 0, 1))
cmip6s$change.cat.maj <- with(cmip6s, ifelse(catmajref == "Arid" & cat.maj %in% c("Cold","Humid"),"Arid to cold or humid",
ifelse(catmajref == "Arid" & cat.maj %in% c("Semi-Arid","Dry subhumid"), "Arid to semi-arid or dry subhumid",
ifelse(catmajref == "Arid" & cat.maj == "Hyperarid", "Arid to hyperarid",
ifelse(catmajref == "Semi-arid" & cat.maj == "Dry subhumid", "Semi-arid to dry subhumid",
ifelse(catmajref == "Semi-arid" & cat.maj == "Arid", "Semi-arid to arid",
ifelse(catmajref == "Semi-arid" & cat.maj == "Hyperarid", "Semi-arid to hyperarid",
ifelse(catmajref == "Semi-arid" & cat.maj %in% c("Humid","Cold"), "Semi-arid to humid or cold",
ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "Dry subhumid to arid, semi-arid or hyperarid",
ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Humid","Cold"), "Dry subhumid to humid or cold",
ifelse(catmajref == "Humid" & cat.maj == "Cold","Humid to cold",
ifelse(catmajref == "Humid" & cat.maj == "Dry subhumid", "Humid to dry-subhumid",
ifelse(catmajref == "Humid" & cat.maj == "Semi-arid","Humid to semi-arid",
ifelse(catmajref == "Humid" & cat.maj %in% c("Arid", "Hyperarid"), "Humid to arid or hyperarid",
ifelse(catmajref == "Cold" & cat.maj == "Humid", "Cold to humid",
ifelse(catmajref == "Cold" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "Cold to dryland",
ifelse(catmajref == "Hyperarid" & cat.maj == "Arid", "Hyperarid to arid",
ifelse(catmajref == "Hyperarid" & cat.maj %in% c("Semi-arid","Dry-subhumid"),"Hyperarid to semi arid or dry subhumid",
ifelse(catmajref == "Hyperarid" & cat.maj %in% c("Humid", "Cold"), "Hyperarid to humid or cold", NA
)))))))))))))))))))
table(cmip6s$change.cat.maj)
##
## Arid to cold or humid
## 2
## Arid to hyperarid
## 491
## Arid to semi-arid or dry subhumid
## 4
## Cold to dryland
## 3
## Cold to humid
## 3700
## Dry subhumid to arid, semi-arid or hyperarid
## 649
## Dry subhumid to humid or cold
## 201
## Humid to arid or hyperarid
## 19
## Humid to cold
## 237
## Humid to dry-subhumid
## 940
## Humid to semi-arid
## 2148
## Hyperarid to arid
## 301
## Semi-arid to arid
## 1462
## Semi-arid to dry subhumid
## 236
## Semi-arid to humid or cold
## 381
## Semi-arid to hyperarid
## 2
colors_changes <- c("Hyperarid to arid" = "#AA44AA",
"Arid to semi-arid or dry subhumid" = "#A9A0F9",
"Dry subhumid to humid or cold" = "#A6D9F9",
"Semi-arid to dry subhumid" = "#aec7e8",
"Semi-arid to humid or cold" = "#377EB8",
"Humid to cold" = "#08306B",
"Cold to humid" = "#ff9896",
"Humid to dry-subhumid" = "#ffbbA9",
"Humid to semi-arid" = "#ffdd55",
"Dry subhumid to arid, semi-arid or hyperarid" = "#CC7733",
"Semi-arid to arid" = "#ff7f0e",
"Arid to hyperarid" = "#d62728")
ggplot() + geom_raster(data = subset(cmip6s, change.catmaj == 1), aes(x=lon, y = lat, fill = change.cat.maj))+
borders()+
scale_fill_manual(values = colors_changes)+
labs(title = "dryer", fill = "")+
theme_minimal()+ylim(-55,90)+
theme(legend.position = "bottom")
5.2.2.1 Changes towards dryer categories
5.2.2.1.1 Historical
map_list <- list()
for(i in c("1850_1880","1985_2015")){
g <- ggplot() +
geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & dryer.maj == "dryer" & !change.cat.maj %in% c("Arid to cold or humid","Cold to dryland","Semi-arid to hyperarid","Arid to semi-arid or dry subhumid")),
aes(x=lon, y = lat, fill = change.cat.maj))+
borders()+
scale_fill_manual(values = colors_changes, na.translate = F)+
labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom", guide_legend(nrow = 2))
map_list[[i]] <- g
}
ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")
5.2.2.1.2 Future
map_list <- list()
for(i in c("2030_2060","2070_2100")){
for (j in c("SSP245","SSP370","SSP585")){
index <- paste(i, j, sep = " , ")
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == j& change.catmaj == 1 & dryer.maj == "dryer" &
!change.cat.maj %in% c("Arid to cold or humid","Cold to dryland","Semi-arid to hyperarid","Arid to semi-arid or dry subhumid","Humid to arid or hyperarid")),
aes(x=lon, y = lat, fill = change.cat.maj))+
# geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & nb.maj > 10), aes(x= lon, y = lat, fill = "'"))+
borders()+
scale_fill_manual(values = colors_changes)+
# scale_fill_viridis_d()+
labs(fill = "Change of category\ncompared to\n1970-2000")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom", guide_legend(nrow = 2))
map_list[[index]] <- g
}}
ggarrange(plotlist = map_list,ncol = 3, nrow = 2, common.legend = T, legend = "bottom", labels = c("SSP 2-4.5", "SSP 3-7.0", "SSP 5.85"))
5.2.2.2 Changes towards wetter categories
5.2.2.2.1 Historical
map_list <- list()
for(i in c("1850_1880","1985_2015")){
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & dryer.maj == "wetter" &
!change.cat.maj %in% c("Arid to cold or humid","Cold to dryland","Semi-arid to hyperarid","Arid to semi-arid or dry subhumid")),
aes(x=lon, y = lat, fill = change.cat.maj))+
borders()+
scale_fill_manual(values = colors_changes, na.translate = F)+
labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom", guide_legend(nrow = 2))
map_list[[i]] <- g
}
ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")
5.2.2.2.2 Future
map_list <- list()
for(i in c("2030_2060","2070_2100")){
for (j in c("SSP245","SSP370","SSP585")){
index <- paste(i, j, sep = " , ")
g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == j& change.catmaj == 1 & dryer.maj == "wetter" &
!change.cat.maj %in% c("Arid to cold or humid","Cold to dryland","Semi-arid to hyperarid","Arid to semi-arid or dry subhumid","Humid to arid or hyperarid")),
aes(x=lon, y = lat, fill = change.cat.maj))+
# geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & nb.maj > 10), aes(x= lon, y = lat, fill = "'"))+
borders()+
scale_fill_manual(values = colors_changes)+
# scale_fill_viridis_d()+
labs(fill = "Change of category\ncompared to\n1970-2000")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom", guide_legend(nrow = 2))
map_list[[index]] <- g
}}
ggpubr::ggarrange(plotlist = map_list,ncol = 3, nrow = 2, common.legend = T, legend = "bottom", labels = c("SSP 2-4.5", "SSP 3-7.0", "SSP 5.85"))
6 Tables
6.1 Changes of aridity category using cat mean
tab <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.AI)[1], Hyperarid = table(cat.AI)[2], Arid = table(cat.AI)[3], Semiarid = table(cat.AI)[4], Drysubhumid = table(cat.AI)[5], Humid = table(cat.AI)[6], total = n())
write.table(tab, "tab.catAI.txt")
tab.future <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(period, model, cat.AI) %>%
summarise(count = n()) %>%
ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
reshape2::dcast(period + model ~cat.AI) %>%
mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
kable(tab.future[order(tab.future$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) | model | period | Hyperarid | Arid | Semi-arid | Dry subhumid | Sum drylands | Humid | Cold | NA | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | historical | 1850_1880 | 3.4 | 6.0 | 7.1 | 3.3 | 19.8 | 28.8 | 48.9 | 2.6 |
| 2 | historical | 1970_2000 | 3.3 | 6.0 | 7.0 | 3.3 | 19.6 | 27.8 | 49.9 | 2.6 |
| 3 | historical | 1985_2015 | 3.3 | 6.1 | 6.9 | 3.3 | 19.6 | 27.3 | 50.6 | 2.6 |
| 4 | SSP245 | 2030_2060 | 3.5 | 6.3 | 7.8 | 3.3 | 20.9 | 28.8 | 47.7 | 2.6 |
| 7 | SSP245 | 2070_2100 | 3.6 | 6.4 | 7.8 | 3.5 | 21.3 | 29.4 | 46.6 | 2.6 |
| 5 | SSP370 | 2030_2060 | 3.2 | 6.7 | 7.6 | 3.0 | 20.5 | 28.8 | 48.1 | 2.6 |
| 8 | SSP370 | 2070_2100 | 4.0 | 6.3 | 8.3 | 3.7 | 22.3 | 28.5 | 46.6 | 2.6 |
| 6 | SSP585 | 2030_2060 | 3.6 | 6.2 | 7.9 | 3.5 | 21.2 | 29.3 | 47.0 | 2.6 |
| 9 | SSP585 | 2070_2100 | 3.8 | 6.4 | 8.4 | 3.6 | 22.2 | 30.1 | 45.1 | 2.6 |
# ggplot(tab.future, aes(x = period))+
# geom_point(aes(y = Hyperarid, col = "Hyperarid", shape = model))+
# geom_line(aes(y=Hyperarid, col = "Hyperarid", group = model, lty = model))+
# geom_point(aes(y = Arid, col = "Arid", shape = model))+
# geom_line(aes(y=Arid, col = "Arid", group = model, lty = model))+
# geom_point(aes(y = get('Semi-arid'), col = "Semi-arid", shape = model))+
# geom_line(aes(y=get('Semi-arid'), col = "Semi-arid", group = model, lty = model))+
# geom_point(aes(y = get('Dry subhumid'), col = "Dry subhumid", shape = model))+
# geom_line(aes(y=get('Dry subhumid'), col = "Dry subhumid", group = model, lty = model))+
# geom_point(aes(y = Humid, col = "Humid", shape = model))+
# geom_line(aes(y= Humid, col = "Humid", group = model, lty = model))+
# scale_color_manual(values = col.cat)+
# theme_minimal()6.2 Changes of aridity category using cat maj
tab.maj <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.maj)[1], Hyperarid = table(cat.maj)[2], Arid = table(cat.maj)[3], Semiarid = table(cat.maj)[4], Drysubhumid = table(cat.maj)[5], Humid = table(cat.maj)[6], total = n())
write.table(tab.maj, "tab.catmaj.txt")
tab.future.maj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(period, model, cat.maj) %>%
summarise(count = n()) %>%
ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
reshape2::dcast(period + model ~cat.maj) %>%
mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
kable(tab.future.maj[order(tab.future.maj$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) | model | period | Hyperarid | Arid | Semi-arid | Dry subhumid | Sum drylands | Humid | Cold | NA | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | historical | 1850_1880 | 4.2 | 7.1 | 8.7 | 0.8 | 20.8 | 28.3 | 48.3 | 2.6 |
| 2 | historical | 1970_2000 | 4.1 | 7.1 | 8.6 | 0.7 | 20.5 | 28.1 | 48.9 | 2.6 |
| 3 | historical | 1985_2015 | 4.0 | 7.0 | 8.4 | 0.9 | 20.3 | 27.3 | 49.9 | 2.6 |
| 4 | SSP245 | 2030_2060 | 4.0 | 7.8 | 9.3 | 0.8 | 21.9 | 28.4 | 47.2 | 2.6 |
| 7 | SSP245 | 2070_2100 | 4.2 | 7.7 | 9.4 | 0.9 | 22.2 | 29.1 | 46.1 | 2.6 |
| 5 | SSP370 | 2030_2060 | 4.1 | 7.2 | 9.4 | 0.8 | 21.5 | 28.6 | 47.3 | 2.6 |
| 8 | SSP370 | 2070_2100 | 4.5 | 7.7 | 10.0 | 0.9 | 23.1 | 28.4 | 45.9 | 2.6 |
| 6 | SSP585 | 2030_2060 | 4.2 | 7.5 | 9.6 | 1.0 | 22.3 | 28.8 | 46.3 | 2.6 |
| 9 | SSP585 | 2070_2100 | 4.5 | 7.9 | 9.4 | 0.9 | 22.7 | 30.0 | 44.6 | 2.6 |
# ggplot(tab.future.maj, aes(x = period))+
# geom_point(aes(y = Hyperarid, col = "Hyperarid", shape = model))+
# geom_line(aes(y=Hyperarid, col = "Hyperarid", group = model, lty = model))+
# geom_point(aes(y = Arid, col = "Arid", shape = model))+
# geom_line(aes(y=Arid, col = "Arid", group = model, lty = model))+
# geom_point(aes(y = get('Semi-arid'), col = "Semi-arid", shape = model))+
# geom_line(aes(y=get('Semi-arid'), col = "Semi-arid", group = model, lty = model))+
# geom_point(aes(y = get('Dry subhumid'), col = "Dry subhumid", shape = model))+
# geom_line(aes(y=get('Dry subhumid'), col = "Dry subhumid", group = model, lty = model))+
# geom_point(aes(y = Humid, col = "Humid", shape = model))+
# geom_line(aes(y= Humid, col = "Humid", group = model, lty = model))+
# scale_color_manual(values = col.cat)+
# theme_minimal()tab.flow <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(period, model, cat.maj) %>%
summarise(count = n()) %>%
ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)
df245 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP245")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>%
subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))
b245 <- ggplot(data = df245, aes(x = period, y = percent))+
geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
geom_label(aes(y = lab.y, label = percent))+
scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
labs(title = "SSP 2-4.5", y = "%", x = "")+
theme_minimal()
df370 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>%
subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))
b370 <- ggplot(data = df370, aes(x = period, y = percent))+
geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
geom_label(aes(y = lab.y, label = percent))+
scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
labs(title = "SSP 3-7.0", y = "%", x = "")+
theme_minimal()
df585 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>%
subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))
b585 <- ggplot(data = df585, aes(x = period, y = percent))+
geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
geom_label(aes(y = lab.y, label = percent))+
scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
labs(title = "SSP 5-8.5", y = "%", x = "")+
theme_minimal()
ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)
6.3 Temporal evolution of percentages, cat.maj
ggarrange(
ggplot(subset(tab.future.maj, period != "1985_2015"))+
geom_point(aes(x = period, y = Hyperarid, col = "Hyperarid", shape = model), size = 3)+
geom_line(aes(x = period, y = Hyperarid, col = "Hyperarid", group = model, lty = model))+
scale_shape_manual(values = c("h","2","3","5"))+
scale_color_manual(values = col.cat)+
guides(col = "none", lty = "none")+
labs(x = "", col = "", shape = "")+
theme_minimal()+
theme(axis.text.x = element_blank()),
ggplot(subset(tab.future.maj, period != "1985_2015"))+
geom_point(aes(x = period, y = Arid, col = "Arid", shape = model), size = 3)+
geom_line(aes(x = period, y = Arid, col = "Arid", group = model, lty = model))+
scale_shape_manual(values = c("h","2","3","5"))+
scale_color_manual(values = col.cat)+
labs(x = "", col = "", shape = "")+
theme_minimal()+
theme(axis.text.x = element_blank()),
ggplot(subset(tab.future.maj, period != "1985_2015"))+
geom_point(aes(x = period, y = get("Semi-arid"), col = "Semi-arid", shape = model), size = 3)+
geom_line(aes(x = period, y = get("Semi-arid"), col = "Semi-arid", group = model, lty = model))+
scale_shape_manual(values = c("h","2","3","5"))+
scale_color_manual(values = col.cat)+
labs(x = "", col = "", shape = "", y = "Semi-arid")+
theme_minimal()+
theme(axis.text.x = element_blank()),
ggplot(subset(tab.future.maj, period != "1985_2015"))+
geom_point(aes(x = period, y = get("Dry subhumid"), col = "Dry subhumid", shape = model), size = 3)+
geom_line(aes(x = period, y = get("Dry subhumid"), col = "Dry subhumid", group = model, lty = model))+
scale_shape_manual(values = c("h","2","3","5"))+
scale_color_manual(values = col.cat)+
labs(x = "", col = "", shape = "", y = "Dry subhumid")+
theme_minimal()+
theme(axis.text.x = element_blank()),
ggplot(subset(tab.future.maj, period != "1985_2015"))+
geom_point(aes(x = period, y = Humid, col = "Humid", shape = model), size = 3)+
geom_line(aes(x = period, y = Humid, col = "Humid", group = model, lty = model))+
scale_shape_manual(values = c("h","2","3","5"))+
scale_color_manual(values = col.cat)+
labs(x = "", col = "", shape = "")+
theme_minimal(),
nrow = 5, common.legend = T, legend = "right"
)6.4 Changes by continent
tab.future.cont <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(period, model, Continent, cat.maj) %>%
summarise(count = n()) %>%
ungroup() %>% group_by(period, model, Continent) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)
tab_list <- list()
for(i in unique(tab.future.maj$Continent)){
tab.i <- subset(tab.future.maj, Continent == i)
df <- tab.i %>%
reshape2::dcast(period + model ~ cat.maj)
tab_list[[i]] <- df
drycats <- which(names(df) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
df <- df %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
df$Cold <- ifelse(is.null(df$Cold) == T, 0, df$Cold)
df$Hyperarid <- ifelse(is.null(df$Hyperarid) == T, 0, df$Hyperarid)
df <- df %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
tab_list[[i]] <- df
}
for(i in names(tab_list)){
df <- tab_list[[i]]
k <- kable(df[order(df$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>%
column_spec(c(8,11), italic = T, include_thead = T) %>% row_spec(2, bold = T)
print(k)
}
for(i in names(tab_list)){
df <- tab_list[[i]]
g <- ggplot(data = df, aes(x=period,y= get("Sum drylands"), col = model))+geom_point()+
geom_line(aes(group = model))+
labs(x="",y = "Percent of drylands", title = i)
print(g)
}
tab.flow <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(period, model, Continent, cat.maj) %>%
summarise(count = n()) %>%
ungroup() %>% group_by(period, model, Continent) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)
df245 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP245")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>%
subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))
b245 <- ggplot(data = df245, aes(x = period, y = percent))+
geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
geom_label(aes(y = lab.y, label = percent))+
scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
facet_grid(rows = vars(Continent), switch = "y")+
scale_y_continuous(position = "right")+
labs(title = "SSP 2-4.5", y = "%", x = "")+
theme_minimal()
df370 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>%
subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))
b370 <- ggplot(data = df370, aes(x = period, y = percent))+
geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
geom_label(aes(y = lab.y, label = percent))+
scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
facet_grid(rows = vars(Continent), switch = "y")+
scale_y_continuous(position = "right")+ labs(title = "SSP 3-7.0", y = "%", x = "")+
theme_minimal()
df585 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>%
subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))
b585 <- ggplot(data = df585, aes(x = period, y = percent))+
geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
geom_label(aes(y = lab.y, label = percent))+
scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
facet_grid(rows = vars(Continent), switch = "y")+
scale_y_continuous(position = "right")+ labs(title = "SSP 5-8.5", y = "%", x = "")+
theme_minimal()
ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)